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Abstract 

A quadrature method for second-order, curved triangular elements in 
the Boundary Element Method (BEM) is presented, based on a polar co- 
ordinate transformation, combined with elementary geometric operations. 
The numerical performance of the method is presented using results from 
solution of the Laplace equation on a cat's eye geometry which show an 
error of order P~^'^, where P is the number of elements. 

1 INTRODUCTION 

The Boundary Element Method (BEM), often called the 'panel method' in fluid 
dynamics, is a standard technique for the solution of boundary integral equations 
in a number of fields. Historically, relatively low-order discretizations have been 
used with geometries modelled using first order elements, and surface variables 
modelled to zero or first order on those elements. There are many methods 
available for the computation of potential integrals on linear panels [l}|6j for 
example], and their behaviour is reasonably well-understood, allowing them to 
be implemented with confidence in production codes. 

More recently, however, there has been increasing interest in the use of 
higher order methods, in part to achieve better geometric fidelity, and in part 
to improve the modelling of solutions. For example, a recently developed panel 
method for whole aircraft aerodynamics |7]|8| , employing accelerated integration 
and summation techniques, depends on the availability of a robust integration 
scheme for second order panels, developed by the authors [9], to avoid some of 
the deficiencies inherent in other curved panel integration techniques (10| for 
example] . 

To clarify the application, we consider the solution of a boundary integral 
formulation of the Laplace equation: 



where 4> denotes potential, x field point position, y position on the surface S, 
and n the outward pointing normal to the surface. The Green's function G is: 

G(x;y) = ^, (2) 
i?=|x-y|. 

To solve this problem using a BEM, the surface S is discretized into a number 
of elements, triangular in this case, over which 4> is approximated by some 
interpolant. This results in a linear system: 

c^W'A.^E/ ^G(x,y)d5,-f / ^^<^(y)d5„ (3) 

where P is the number of elements (panels), i is the index of a surface point, 5^- 
is the surface of panel j, and the constant Ci(x) is a geometric property given 
by: 

, , f gG(x^,y) , 

equal to 1/2 at a smooth point on the surface, taking some other value at sharp 
edges. Inserting the interpolant for each element into Equation [3] yields a system 
of equations relating (j) to d(j)/dn, allowing the problem to be solved subject to 
the specification of some boundary condition. In aerodynamic problems, this 
will usually be the Neumann boundary condition, specifying the surface normal 
velocity d4>/dn. Upon solving for surface potential (j), the boundary integral 
can then be used to compute the potential or its derivatives, i.e. fluid velocity, 
external to the surface, using Equation [l] 

The core of the implementation is then the evaluation of the panel integrals 
in Equation [3) For planar panels, there is no great difficulty, and for the Laplace 
equation, the integration can be performed analytically using a variety of ap- 
proaches [l}|6] , although a numerical method is still necessary for the Helmholtz 
equation for acoustic scattering and radiation. When the panel is curved, how- 
ever, a fully numerical method is required, and extra difficulties arise in finding 
a transformation which maps the curved panel to a reference domain where 
standard quadrature rules can be applied. 

A common approach in integrating over planar panels is to convert to po- 
lar coordinates with axis perpendicular to the element plane, which mitigates 
diflticulties caused by the 1/R singularity in the Green's function. Such an ap- 
proach has been used in computing the self-term for curved panels [9], and a 
similar technique will be used here. Alternatives which have been used include 



the mapping of the element onto a plane triangle 10 , which can, however, only 
be used for the single layer potential, and onto a sphere [9], with appropriate 
conversions between the appropriate Jacobians. 

In this paper, we present a technique for the evaluation of a quadrature rule 
for second order curved panels which uses a polar transformation of the integral 
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Figure 1: Description of second-order curved triangle. Left hand image: curved 
triangle; right hand image: reference right-angle triangle. 

combined with basic geometric operations, to give a method more akin to the 
current techniques for planar panels, but with additional complexity due to the 
need to perform the geometric operations on curved element edges. 

2 ANALYSIS 

The quadrature method for the curved element is derived for a panel in a refer- 
ence position. For a planar element, there is no difficulty in defining an element 
plane which can be used to fix a coordinate system for the triangle, but for 
the curved elements we consider here, there is clearly a choice to be made. A 
standard approach is to use a plane tangent to the element, through some appro- 
priate point, for example, the field point x when a self-term is being computed, 
but here we use the plane defined by the three corners of the triangular element. 
The problem axes are rotated and shifted so that the corners lie in a plane z = 
and the field point x = (0, 0, z), i.e. the origin of the coordinate system is taken 
as the projection of the field point onto the triangle reference plane. In this ori- 
entation, cylindrical polar coordinates can be readily defined and used to carry 
out the required integration. The rest of this section describes the geometric 
operations employed, and the technique used to define the quadrature rule. 

2.1 Description of second order triangles 

Figure [l] shows the notation used for description of the second order triangle. 
Nodes are numbered 1,2,3 for the corners, and 4,5,6 for the points internal to 
an edge. Quantities on the element, including position, are interpolated using 
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the second order shape functions for the reference element: 

'.U{£.,il)yi, (5) 
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0<C<1, 0<ry<l-C- 

As shown in Figure [T] the edges of the triangle are defined by second order 
interpolation on three points, and can be described using a single variable 7, < 
7 < 1, with 7 increasing in the anti-clockwise direction on each edge. Inserting 
the conditions for each edge shown in Figure [l] gives three shape functions: 

Ji ^ 272 -37 + 1, (8a) 

J2 = 272 - 7, (8b) 

•h = -47' + 47, (8c) 

with the edge described by: 

y(7) = ytJiin) + YjMi) + Vi+sMi), (9) 

where (i, j)=(l, 2), (2, 3), (3, 1) for each edge respectively. A point on an edge 
can be given in the general coordinates (^, rj) using the relations: 

r (7,0) t = l; 

i^,v) = { (1-7,7), *-2; (10) 

[ (0,1-7), « = 3. 

These shape functions will prove useful in determining intersections between 
edges and lines in the plane, necessary in finding the domain of integration for 
quadrature over the triangular element. 

2.2 Integration over the triangle 

Figure [2] shows the curved triangle in its reference position, rotated so that 
the corners lie in the plane z — 0, and shifted so that the field point lies at 



4 



z 




Figure 2: Coordinate system for calculation of integrals. The triangle is rotated 
so that the vertices 1,2,3 lie in the plane z = (shown shaded), and translated 
so that the field point lies at a position (0, 0, z). 



X = (0, 0, z). The integral to be evaluated is 

t t ^ f{tv)Jitv)dvd^, (11) 
Jo Jo 

where J(^, rj) is the Jacobian for the transformation from (^, t]) to coordinates 
on the element surface. Converting to Cartesian coordinates in the problem 
system of axes: 

1= f f{^,v)rdrd9, (12) 

J A 

where integration takes place over the projection of the curved triangle into the 
plane z = and the function /(C^y) is computed by transformation from {r,6) 
to (e,?7). 



The integration of Equation [12] is conceptually simple, and is readily applied 
to first order elements [6], but gives rise to some extra complexities in the 
second order case, shown in Figure [3| As written, the integration is composed 
of a sequence of integrals over r, along rays at fixed values of 0. In Figure |3j two 
such rays are shown, at angles 9i and 62- The ray at 9i presents no particular 
difficulties: it has one entry and one exit point on the element boundary, both 



easily found using analytical methods (see Section 2.3). The ray at 02, however, 
is broken as it crosses the triangle boundary, having two entry and two exit 
points. In evaluating Equation |12[ this case must be handled, as must the case 
of a ray which lies tangent to an edge. The algorithm of Section |2.4| handles 
these special cases, using the geometrical operations of the next section. 

2.3 Geometrical operations 



The quadrature algorithm, presented in Section 2.4 depends on the availability 



of a number of elementary geometrical operations, described in this section. 
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Figure 3: Intersection of rays with curved triangle: bold lines show the part of 
each ray which lies inside the domain of integration. 

These operations can be implemented analytically using standard methods and 
are used to determine the limits of integration in 6, to break the integral at 
possible points of discontinuity, and in r, to find the entry and exit points on 
the triangle boundary. 

The first operation is finding the intersection, 7 or (r, 6*), between a ray 
through the origin of angle 9 and an edge of the triangle. The coordinates of a 
point on the edge are given by Equation [9} 

X = XtJii'j) + XjJ2{i) + Xi+3 J3(7), y = ViJiij) + yjJ2{j) + Vi+sJail), (13) 

while the ray is given hy x — r cos 6, y — r sin 6, so that, upon substitution: 

{xi sin 9 — yi cos 9)Ji{'-f) + (xj sin 6* — yj cos (?) 72(7) + (a^i+s sin 9 — yi+y, cos9)J3{j) = 0, 

(14) 

which is a quadratic in 7. To find the intersection: 

1. solve Equation [T4| for 7; 

2. for each value of 7, with < 7 < 1 

(a) compute (x, y) and r = xj cos 9 ox r = yj sin 9\ 

(b) if r > 0, 7 is a valid intersection. 

In this operation, r is computed as shown in order to accept only r > 0, to avoid 
double counting of intersections. The condition < 7 < 1 is imposed to exclude 
corners of the triangle, as these are handled separately in the algorithm. 

Tangents which pass through the origin must also be determined, in order 
to break the integration at these points. They are found using the equation of 
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a tangent to a point {xo,yo) on the edge: 



y = yo + 



dy 

dx 



{x - Xo), 



dy 

dx 



yiJ'ii.i) + y]J'2ii) + Vi+aJ'sii) 



where the prime denotes differentiation with respect to 7. Setting x = y 
find a tangent through the origin yields: 



(15) 
(16) 
to 



Xiy, [J'i{l)J2{l) - J2{l)Ji{l)\ + a^i2/»+3 [^2(7)^3(7) - J'MUl)] 

+ Xi+zyi [4(7) (7) - (7)^3(7)] = 0, (17) 

which is a cubic which can be solved for 7 subject to the constraint < 7 < 1 
and that 7 be real. 

The final part of determining the limits of integration in 8 is the angle t/j of 
a tangent to an edge, given by: 



ijj = tan" 



-1 yiJijj) + VjJfiii) + y»+3 ^3(7) 

XiJ[{l) + XjJ'^i^) + Xj+s J^(7) ' 



(18) 



from the slope of the curve at a point 7 on an edge. A second angle V' + tt is also 
included in order to ensure that rays in both tangent directions are included. 



2.4 Quadrature algorithm 

The algorithm for the quadrature rule consists of a first stage in which the range 
of integration in 9 is broken into a set of intervals, and a second in which the 
integration is performed over these intervals. Initially it must be determined 
whether the origin lies inside, outside or on the boundary of the triangle pro- 
jected into the plane z = Q. This is done by first checking if the origin lies 
outside a box containing the points (2:2,2/2), and {x^,yQ). If not, its 

coordinates rf) in the triangle are found using Newton's method, and it is 
checked whether they lie within the triangle. 
In the first stage of processing: 

1. find possible limits of integration 0^, i = 1, . . . , A'', as the angles of rays 
joining the origin to the triangle's corners, the angles of tangents through 
the origin, and, if the origin lies on an edge, the angles of tangents to the 
edge; 

2. adjust all angles to lie in the range 0, 27r; 

3. sort the list of limits Oi in ascending order; 

4. if the origin lies inside the triangle, append the angle 0i + 2tt. 

Given the list of N angles from the first stage, the nodes and weights of the 
quadrature rule are found for each pair of limits, 9i, 6i+i by this procedure: 



7 



1. select a quadrature rule with abscissae 9k and weights w^, k — 1, . . . , ii'; 

2. for each k: 

(a) find the radii Vj of the intersections of a ray of angle Ok with the 
triangle edges (prepend r = if the origin lies inside the triangle); 

(b) for each pair of limits r2j and ^■2^+1, select a quadrature rule (r„, wJJ^), 
m — 1, . . . , AI; 

(c) for each point (r^ cos 9k, rm sin 6'^): 

i. find the corresponding coordinates 77) using Newton's method; 

ii. append to the quadrature rule the abscissa rj) and the weight 
fmW^w^j^/ J^i^^Trf), where J2 is the Jacobian for conversion from 
(C,r?) to {t,9). 

The resulting quadrature rule can be used to evaluate an integral on the panel 
by summation: 

/(^, V)J{^, V) dvd^-Y. Vn)J{in, nn)Wn, (19) 

n 

where J{^,r]) is the Jacobian for conversion from (^,77) to the surface coordi- 
nates (x,y,z), allowing the rule to be used in the same manner as standard 
quadratures for triangles. We note that since the quadrature rule is mapped to 
the reference triangle, the sum of the weights Wn should be equal to 1/2, which 
gives a convenient error measure for checking the accuracy of the quadrature. 

Finally, a good starting guess for Newton's method in finding coordinates on 
the reference triangle is provided by the intersection point (r, 9) since its value 



of 7 on the edge is known, and can be converted to 77) using Equation 10 
Each set of coordinates (^, ry) on a ray can then be used as an initial guess for 
the evaluation at the next quadrature point. 



2.5 Quadrature selection 

The quadrature rule is implemented using a sequence of Gaussian quadratures 
for integration in 9 and r. This has the first advantage that the endpoints of the 
integral are not included. Since tangents are used to fix limits of integration, 
there is then no ambiguity in determining the number of entry and exit points 
in radius. 

The quadrature rules are selected using a criterion which gives some adaptiv- 
ity to the intervals of integration. A point separation parameter A9 is specified 
and used to estimate the number of points required in the quadrature rule: 

K = - 9,)/A9] , (20) 

where [ • ] denotes rounding of the value to the nearest natural number. The 
resulting value K is then adjusted to lie in a range ifmin ^ K < K^i^-^. A similar 
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method is used to select quadrature rules in r. In the calculations presented 
here, the values of A6' and Ar are computed with user-defined constants Ne and 
Nr, and set as follows: 

A^=^, (21a) 

Ar=^^, (21b) 

where (pi is the angle subtended by the corner of the triangle at vertex i, i = 
1, 2, 3 and £i is the length of the straight edge starting at corner i. This gives a 
quickly computed abscissa density under user control. 

Finally, in using the algorithm in a code, a criterion is required in deciding 
when to use it, and when not. In this case, a parameter a is computed based 
on easily evaluated geometric properties of the element. These are the mean 
values of the nodes and the radius of a sphere containing the element: 



y = ^Ey- (22a) 

p = smax|yi-y|, (22b) 

i 

where s is a scaling factor, with, in this case, s = 2^/^. Given these values for 
the element, a is determined as follows: 

1. if X lies on the element, cr = 0; 

2. compute = |x - y|: ii Px > p, cr = p^/p; 

3. otherwise, a = z/p, 

where z is the distance of x from the element reference plane, as noted above. 
This gives a quickly-computed parameter which varies from on the element 
to large values away from the element, but remains small in some reasonable 
neighbourhood, so that it can be used to select quadrature methods. 



3 NUMERICAL TESTS 

The algorithm of the previous section is demonstrated using two sets of results. 
The first is an illustration of the distribution of quadrature nodes, using a sample 
element, while the second is an assessment of the accuracy and convergence of 
the method when implemented in a BEM program. 

3.1 Quadrature points 

To demonstrate the nature of the quadrature point distributions generated by 
the algorithm, an element with one straight, one convex and one concave edge 
has been used, with the origin placed inside and outside the element, on a 
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Figure 4: Quadrature points for different field points: origin is slrown as a 
cross, quadrature points as dots. Top row: origin inside element, on vertex, and 
outside element. Bottom row: origin on straight, convex and concave edge. 
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Figure 5: The cat's eye geometry 



corner, and on each of the edges in turn. In order to show the point distribution 
more clearly, fixed length quadrature rules have been used, with sixteen points 
in both angle and radius. This gives a higher than normal density in small 
intervals of integration, and a lower density in larger regions, but is helpful for 
visualization. 

Figure |4] shows the resulting quadratures, with the origin indicated by a 
cross. Each of the cases considered gives rise to a qualitatively different point 
distribution. With the origin located inside the element, the region of integration 
is divided by rays to each of the corners (compare Figure 8 of reference [o]). 
When the origin is placed on a corner of the element, there are two clearly 
demarcated domains of integration, separated by the tangent to the curved 
edge at that corner. Similarly, when the origin is moved outside the element, 
there are two domains, separated in this case by a ray joining the origin to the 
most distant corner. 

When the origin lies on an edge, the situation is slightly more complicated. 
When it is on the straight edge, the element is divided into two regions, separated 
by the ray to the furthest corner. When it lies on the convex edge, there is a 
thin region of integration bounded by the edge and the rays to the vertices on 
that edge. Conversely, on the concave edge, the narrow region of integration is 
bounded by the edge and the tangent to the edge at the origin. 

3.2 Numerical accuracy and convergence 

The accuracy and convergence of the integration method are tested by imple- 
menting it in a BEM code flTl and solving the Laplace equation on a cat's eye 
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geometry, shown in Figure [5] This is a unit sphere with one octant removed 



recommended as a more stringent test of BEM codes than a simple sphere 12 



since it contains discontinuities in the geometry, as would be found, for example, 
in aerodynamic calculations. The surface was meshed using GMSH 13 , chang- 
ing the discretization length to produce panels of varying sizes. A second mesh 
was produced for comparison by splitting each second order element into six 
triangles, giving a mesh of planar elements based on the same nodes. The polar 
quadrature rule was selected for cr < 1, with Ng — Nr = 8, and a twenty-five 
point symmetric rule 14 for a > 1. 

A Neumann boundary condition was generated using a point potential source 
positioned inside the surface at (—0.2, —0.2, —0.2). Solving Equation [l] for the 
surface potential (j) gave a result which could be compared to the result found 
analytically for the point source. The error estimate is the r.m.s. difference 
between the computed and the analytically specified data. The error is plotted 
in Figure [6] as a function of element number, and, for convenience, of node 
number. The error is fitted using a power law, to estimate the convergence rate, 
and the superior numerical performance of the second order method is clear. 
Plotting against panel number shows a similar convergence rate as in other 
work [91 Figure 12], for linear panels, and P~^/^ for quadratic. Plotting 
against node number, which can be taken as a proxy for the memory requirement 
for the matrix used to solve the problem, and which shows error for different 
element types on the same point distribution, shows similar trends. 



4 CONCLUSIONS 

A quadrature technique for second order triangular elements has been presented 
and tested on a realistic geometry. It has been found that the method is accurate 
and convergent, giving an error which scales as P^^-^ in the example tested. It 
is concluded that the technique is readily implemented and can be used as a 
direct replacement for existing quadratures. 
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